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O , Abstract 

o . 

, We present a simple formalism for parton-shower Markov chains. As a first step towards more 

^ ' complete 'uncertainty bands', we incorporate a comprehensive exploration of the ambiguities in- 

1-^ , herent in such calculations. To reduce this uncertainty, we then introduce a matching formalism 

• which allows a generated event sample to simultaneously reproduce any infrared safe distribution 

CNl , calculated at leading or next-to-leading order in perturbation theory, up to sub-leading corrections. 

' To enable a more universal definition of perturbative calculations, we also propose a more general 

(— I I definition of the hadronization cutoff. Finally, we present an implementation of some of these ideas 

pi i' for final-state gluon showers, in a code dubbed ViNClA. 

r^/ 1 Introduction 

^ I At present collider energies, the strong (QCD) coupling strength is sufficiently large that even the 

^ ■ most sophisticated approximations are typically reliable only over a limited region of phase space. De- 

lO . scriptions which work well for "hard" radiation (extra jets, hard bremsstrahlung) break down rapidly in 

soft and/or coUinear regions (jet structure, soft bremsstrahlung), and vice versa. In addition, at scales 
below a GeV or so, non-perturbative effects must also be taken into account, and a transition made 
to a description in terms of screened charges inside colorless hadrons. In this paper, we shall put the 
I main focus on aspects which are systematically examined in two complementary perturbative approx- 

O ■ imations: that of fixed-order truncations, appropriate for "hard" radiation, and that of parton shower 

t> . resummations, appropriate for "soft/coUinear" radiation. 

^ , In the past, these two approximations were often pursued independently. The last decade or so 

?H I has witnessed rapid progress in our understanding of how the virtues of each can be used to overcome 

■ — ' the vices of the other, to yield "matched" results which attempt combine the best features of both ap- 

proaches. The earliest concrete approach, due to Sjostrand and collaborators [1], is implemented in the 
Pythia generator [2] and consists of re-weighting the first parton shower emission off a hard system 
X by a correction factor equal to the ratio needed to reproduce the tree-level {X + 1) -parton matrix 
element. This reweighting relies heavily on the shower algorithm(s) in Pythia [3] covering all of phase 
space and on the first shower emission being clearly identifiable as the "hardest". A different approach 
was needed for the coherent angular-ordered algorithm [4] used in the Herwig Monte Carlo [5], in 
which soft gluons can be emitted at large angles "before" harder ones at smaller angles. The approach 
developed by Seymour for this puipose [6] combines two ingredients: in the region populated by the 
shower, emissions are reweighted to produce the matrix element rate almost as above (generalized to the 
angular-ordered case), whereas in the so-called "dead zone," separate {X + l)-parton events are gener- 
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ated according to the appropriate matrix element, and weighted by a Sudakov factor. This technique can 
be viewed as a precursor to the modern CKKW approach [7]. 

Although matching beyond one extra parton was attempted within the Pythia framework [8], the 
complexity of the problem grows rapidly. The CKKW and MLM [9] matching schemes broke through 
this banier, in principle providing a framework for matching through any number of tree-level matrix 
elements, though practical applications are still limited to including matrix elements for at most a hand- 
ful of additional partons. The original CKKW approach is implemented in the Sherpa generator [10]. 
The MLM one is in principle less dependent on the specific implementation, but is probably most of- 
ten used with Alpgen [11] interfaced [12] to either Herwig or Pythia. The basic idea behind the 
CKKW scheme has since been refined and extended, first via so-called pseudoshowers introduced by 
Lonnblad [13, 14] in the context of the color dipole model [15, 16] and implemented in the Ariadne 
generator [17]; and later by Mrenna and Richardson [18] using MadGraph [19], again interfaced 
to Herwig and Pythia. The most recent advance to be implemented in a widely-used generator 
was the subtraction-based loop-level matching proposed by Frixione and Webber [20] that led to the 
Mc@NLO [21] add-on to the HERWIG generator. 

More recently, several groups have presented proposals to improve Mc@NLO-style matching [22]; 
to include one-loop contributions in a CKKW-like scheme [23]; to develop a formalism capable of 
dealing with subleading color and spin effects [24]; to include small-x effects [25]; and to use Soft- 
CoUinear Effective Theory (SCET, see ref. [27] and references therein) as a framework for matched 
parton showers [26]. 

Making use of a generalization of the matching proposed by Frixione and Webber, our aim is to 
present a simple formalism for leading-log leading-color parton showers, constructed explicitly with 
two main goals in mind: 1) including systematic uncertainty estimates, and 2) combining the virtues of 
CKKW-type matching (matching to tree-level matrix elements with an arbitrary number of additional 
partons) with those of Mc@NLO-type approaches (matching to one-loop matrix elements). This can be 
done in a manner which is simultaneously simple and does not introduce any dependence on clustering 
schemes or p± cutoffs beyond those required to regulate explicit subleading logaiithms and hadroniza- 
tion effects. Negative weights will arise in general, in particular when matching to one-loop matrix 
elements. Within the formalism we present, however, they can typically be avoided through judicious 
choice of subtraction terms. Moreover, phase space generation only needs to be carried out on matrix 
elements which have had their singularities subtracted out, and hence should be relatively fast. Finally, 
as a first step towards making hadronization models (in paiticular their "tuning") less dependent on the 
details of the parton shower they are used with, we propose a generalized definition of the hadronization 
cutoff. 

As an explicit proof of concept, we have combined the antenna factorization formalism [28] with 
that of dipole showers [16] in a code dubbed ViNCiA (Virtual Numerical Collider with Interleaved 
Antennae) [29], which is being developed both stand-alone and as a plug-in to Pythia 

The organization of the paper is as follows: in section |2l we define the general formalism for parton 
showers, including a brief discussion of each component. We then expand the shower into partial cross 
sections of fixed multiplicities of resolved paitons. This expansion is used in section [3] to derive a 
set of matching terms for abitrary tree- and one-loop matrix elements, up to corrections of order an 
infrared cutoff (hadronization scale) and subleading logarithms. Section|4]then deals with improving the 
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infrared factorization between parton showers and hadronization models. The numerical implementation 
of many of our ideas, in the form of the ViNCiA code, is then presented in section [5] In section [6l we 
discuss how one might go further in the perturbative expansion. We then round off with conclusions and 
outlook in section |2l 



2 The Shower Chain 



As a starting point, consider a Markov chain algorithm [30] ordered in some measure of time t. Such 
chains characterize the development of a broad variety of systems. In our application, the system will 
be a set of partons and the role of t will be played by a measure of parton resolution scale Q = 1/t, but 
the chain could equally well represent the real-time evolution of a simple system such as an ensemble 
of radioactive nuclei. Two global quantities characterize such evolution: the starting configuration and 
the duration of the experiment tend — ^in- At each differential time step in between, there is a probability 
density A(t) for the system to undergo a non-trivial change. After such a change, A{t) itself may 
change (for example, a nucleus may be replaced by its decay products). For the chain to have the 
Mai'kov property, all that is required is that A{t) depend only on the system's present state, not on its 
past history. This property will turn out to be useful in the context of higher-order matching in section 

m 

Let us write the Sudakov factor [31], the probability that the system does not change state between 
two times tin and tend, as: 

A(tin, tend) = exp (- [ ' dt A{t)] , (1) 



iin 

with the understanding that A can depend on the particular system configuration at time tin, and thereby 
that A implicitly has such a dependence as well. 

In the parton shower context, A{t) is the total paiton evolution or splitting probability density. This 
includes sums and integrals over all possible types of transitions, such as gluon-to-gluon pair or gluon- 
to-quark pair splitting. As our aim is to resum the leading singularities, A{t) must necessarily be infrared 
divergent: there is an infinite probability to radiate a soft or collinear gluon. The evolution variable must 
therefore itself be infrared safe, such that all the singularities of A correspond to "late times" as defined 
by t (or, equivalently, 1 / Q). Specifically, for leading-log evolution t must regulate at least all leading-log 
divergences. 

Letting {p}n denote a complete specification of an n-parton configuration in the leading-color limit 
(carrying not only information on momenta, but also the color ordering, flavors, and perhaps polariza- 
tions), we define the leading-log Sudakov factor by. 



/ /-tend f d$i^L r-, 

A(tn,tend;Mn) = expl - / dt„+i / — 5(t„+i - 1'*' ) ) 



(2) 

where d<I>„ denotes the n-particle Lorentz-invariant phase-space measure, so that d$|^!|_]^/d<I>„ repre- 
sents the branching phase space, and the t-ordering is imposed via the integral over t„-|_i together with 
the 6 function. This definition will be the cornerstone for the remainder of this paper, and so we now 
devote a few paragraphs to its explanation. To simplify the notation, we shall usually let the dependence 
on {p}n be implicit, letting A(tn, tend) = A(tn, tend; {p}n)- 
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The first important aspect is that we define the Sudakov factor not for a lone parton or color dipole, 
but rather for the n-parton configuration as a whole. The degree to which the evolution of smaller 
subsystems factorize will naturally play an important role but does not need to be specified explicitly 
at this point. The sum over i E {n — > n + 1} runs over all possible ways of obtaining n + 1 partons 
from the original n ones. For example, starting with an n-parton configuration of which Ug partons are 
quarks, an unpolarized parton shower with four active quark flavors would yield one term for each quark 
in the event (n^ x q ^ qg) and five terms for each gluon {n — Uq) x (jif x g ^ qq + g ^ gg). 

The evolution phase space is represented symbolically by the (n + 1) -parton phase space for an 
evolution step of type i, d<I>„i^]^, divided by that of the evolving configuration, d<I)„. Its specification 
requires three variables, along with a mapping from these variables to the phase space for the emission. 
Existing parton-shower Monte Carlo implementations each choose a different function for this map. 
To name a few known issues, the map may have 'dead zones' where it is zero, and/or it may have 
regions where several independent emitters i populate the same point (such double counting is not 
necessarily a problem, so long as the sum is properly normalized); in some formulations, the entire 
event may participate in each branching, in others only a specific pair of partons 'recoil' off each other; 
in analytical leading-log (LL) resummations, a purely collinear map is usually used, which slightly 
violates momentum conservation, and so on. Our point is not which choice is 'best', but that many are 
possible, each leading to a different shower evolution. 

In this paper, we shall require that the partons at each step be on shell and that four-momentum be 
conserved. An (n + 1) -parton phase space then has 3(n + 1) — 4 degrees of freedom: 

/n+1 n \ n+1 ,3_. 

d*„,,=(2.)'*'(g«-gB]n^, (3) 

where we denote the n original momenta by p, and the notation for four- vectors is pi = [Ei, pj). This 
represents all possible (n+l)-paiton configurations consistent with energy and momentum conservation. 
In the context of evolution, however, we are already implicitly integrating over all possible n-parton 
configurations, and we are also explicitly summing over all possible evolution possibilities i for each 
such configuration. The notation d$|^!^]^/d<I>„ is thus intended to signify the subdivision of the full 
phase space into discrete (but possibly overlapping) regions, each corresponding to a specific n-paiton 
configuration and a specific i. Since the n-parton phase space has 3n — 4 degrees of freedom, this 
amounts to imposing 3n — 4 additional 6 functions and adding the explicit superscript [i] : 

^ = (e ft- - Eft) n - {f ).) ■ (4) 

where k^"^ is the inverse of the map discussed above. This inverse can be viewed as a clustering defi- 
nitioiil that, given i and the (n + 1) -parton configuration, reconstructs the corresponding "unevolved" 
one, {p}n- The requirements on the map are thus: 

1. For each i, a unique inverse of ki must exist {ki must be injective), however each individual Kj 
does not necessarily have to cover all of phase space, since we only care about the coverage after 
summing over i and integrating over d<I>„. 

"That is, "inverts" the shower in a manner similar to the action, e.g., of Arclus on the ARIADNE shower [17] or 
Pyclus on the px-ordered Pythia shower. 
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2. After summing over i and integrating over d<I>„, the resulting composite map should cover all of 
phase space (be surjective), in order to avoid creating dead zones. It does not necessarily have 
to be one-to-one: the (n + l)-parton phase space may be covered several times so long as this is 
properly taken into account in the normalization of the radiation functions Ai (or more precisely, 
of their singular parts). 

Obviously, these statements only apply to configurations that are supposed to be obtainable via shower 
branchings in the first place, and not, for instance, to subleading color topologies such as Z — > qqgg 
with the two gluons in a color singlet state. 

Below, we shall restrict our attention to maps corresponding to dipole-antenna showers, such that: 

Dipole Showers: = pl^^ . (5) 

The first factor on the right hand side indicates that we we choose 3(?i — 2) 6 functions to express 
that n — 2 partons don't move at all. This leaves us with a dipole-antenna phase space, d$3Vd$2> 
carrying nine degrees of freedom compensated by the six remaining S functions. Introducing the notation 
a + b^a + r + b which we shall use for dipole-antenna branchings throughout, four of these delta 
functions embody overall momentum conservation, 

S'-^HPa + Pr + Pb - Pd - Pi) ■ (6) 

The last two delta functions specify the global orientation of the plane spanned by the three daughter 
partons in the center of mass of the parent dipole (the branching plane), relative to the axis of the parent 
dipole, in terms of two angles, 9 and ip. Parity conservation fixes one of them so that the branching 
plane contains the parent dipole axis: 

6{e-9) , (7) 

where 9 is the orientation angle of the parent dipole in a global coordinate system. The S function in 
fixes the rotation angle of the daughters around an axis perpendicular to the branching plane. This 
angle does have a reparametrization ambiguity away from the collineai^ and soft limits and hence has the 
following general form 

S{iP - iP - iPad) , (8) 

where is the other global orientation angle and the reparametrization term tpad will be explored further 
in section |5] (see also ref. [32]). For the time being, we note only that it must vanish in the soft and 
collinear limits. 

The remaining three integration variables we shall map to two invariant masses, s^r = {Pa + PrY 
and Srb = {Pr +Pb)^> and the last Euler angle, (p, describing rotations of the branching plane around the 
parent dipole axis. The antenna phase space then takes the following form: 



1 



A(sW,m?,m?) ' dcf> diljdcos9 

— — dSnrdSrh (") 

d$2 167r2 "''27r d^/'dcos^ 

= \ m (^Sardsrb^ for mi, = mT = ; ip = ip + ipad and 9 = 9, (10) 
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where A(a, b, c) = + 6^ + — 2ab — 2bc — 2ca is the Kallen function, s'*! is the invariant mass squared 
of the branching dipole, and m. g are the rest masses of the original endpoint partons. The second line 
represents the massless case, with the two orientation angles 6 and ip fixed as discussed above. 

Immediately following the phase space in eq. (O is a 5 function requiring that the integration variable 
tn+i should be equal to the ordering variable t evaluated on the set of n + 1 partons, {p}n+i, i.e. that the 
configuration after branching indeed corresponds to a resolution scale of tn+i- We leave the possibility 
open that different mappings will be associated with different functional forms for the post-branching 
resolution scale, and retain a superscript on t'*! to denote this. 

Finally, there are the evolution or showering kernels Ai{{p}n—>- {p}n+i), representing the differen- 
tial probability of branching, which we take to have the following form, 

^i(Mn^Mn+i) = 47ras(/ii?({p}„+i)) Ci ai({p}„^{p}„+i) , (11) 

where Airas = gl is the strong coupling evaluated at a renormalization scale defined by the function 
^R, Ci is the color factor (e.g. Ci = Nc = 3 for gg —>■ ggg), and is a radiation function, giving a 
leading-logarithmic approximation to the corresponding squai^ed evolution amplitude (that is, a paiton 
or dipole-antenna splitting kernel). When summed over possible overlapping phase-space regions, the 
combined result should contain exactly the correct leading soft and coUinear logarithms with no over- or 
under-counting. Non-logarithmic ('finite') terms are in constrast arbitrary. They con^espond to moving 
around inside the leading-logarithmic uncertainty envelope. The renormalization scale fiji could in 
principle be a constant (fixed coupling) or running. Again, the point here is not to impose a specific 
choice but just to ensure that the language is sufficiently general to explore the ambiguity. 

Together, eqs. (HJl, and (fTTI) can be used as a framework for defining more concrete parton 
showers. An explicit evolution algorithm (whether based on partons, dipoles, or other objects) must 
specify: 

1. The choice of perturbative evolution variable(s) tt'l 

2. The choice of phase-space mapping d$[*!^]^/d$„. 

3. The choice of radiation functions Oj, as a function of the phase-space variables. 

4. The choice of renormalization scale function 

5. Choices of starting and ending scales. 

The definitions above are already sufficient to describe how such an algorithm can be matched to 
fixed order perturbation theory. We shall later present several explicit implementations of these ideas, in 
the form of the ViNCiA code, see section [5] 

Let us begin by seeing what contributions the pure parton shower gives at each order in perturbation 
theory. Since A is the probability of no branching between two scales, 1 — A is the integrated branching 
probability ^branch. Its rate of change gives the instantaneous branching probability over a differential 
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time step dt„+i: 



din+1 dt„+i 



1 — 



(12) 

This expression still contains an explicit integral over all phase-space variables except tn+i- The corre- 
sponding fully differential distribution of the (time-ordered) branching probability is simply the "naive" 
evolution kernel times the Sudakov factor, 

Ai{{p}n^{p}n+l)^{tn,tn+l) . (13) 

We now seek the complete fixed-order expansion for the distribution of an observable O, as com- 
puted by the Markov process. By definition, O is always evaluated on the final configuration, reached 
once the Markov chain terminates at tend- A convenient physical interpretation of tend is as the hadroniza- 
tion cutoff, thad> beyond which the parton shower is not evolved; since the evolution would receive 0{l) 
coiTcctions from hadronization beyond this scale, exclusive event properties can only be further probed 
by switching to a non-perturbative description in that region (alternatively, stopping the shower and 
evaluating O on the partons at this scale is still an improvement over fixed-order calculations). We will 
return to a generalization of this notion in a later section; for now merely take thad as a cutoff in the 
evolution variable. This amounts to implicitly treating all radiation (and hadronization effects) below 
the cutoff inclusively, that is summing over additional partons below that scale. 

Starting from any process, X, with differential phase space weight wx, the parton-shower improved 
distribution of O is: 

^|pg = / d'^xwxS{{p}x,0) . (14) 

The definition of wx in the context of matching will be explored in the next section, but for the pure 
shower wx is just the tree-level matrix element squared for the parent process X, possibly subject to 
matrix-element-level cuts or constraints (e.g. Z production restricted to a window around the Z mass, 
etc). S" is a showering operator that generates the Markov chain starting from a list of paitons {p}x- It 
is defined by: 

S{{p}x,0) = 5(0-0({p}x))A(tx,thad) 



X + exclusive above 1 / thad 
+ 'dtx+iEy -^^S{tx+i-t^'^{{p}x+i))A{tx,tx+i)A{...)S{{p}x+uO) 



X + I inclusive above 1 / thad 

(15) 

where tx, the starting scale for each successive step of the evolution, depends implicitly on {p}x, the 
integration over tx+i runs over all possible branching scales between tx and thad> and Ai{...) is defined 
by eq. ([TTI ). Expanding the Markov chain to a few orders will be useful in the context of matching below 
and simultaneously illustrates expUcitly how the chain works: 
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sxmx,o) 



5(0-0({p}x))A(tx,thad) 

., ' 



X + exclusive above 1 / thad 




6{tx+i - t^'^{{p}x+i))A{tx,tx+i)A{...) 



X A{tx+1, thad)5 (O - 0{{p}x+i)) 



(16) 




X + 1 exclusive above 1 / thad 



X \ A{tx+2, thad)'^ {O - 0{{p}x+2)) + - ■ ■ 



X + 2 exclusive above 1 / thad 



Each underbraced term corresponds to the finite contribution from a specific exclusive final-state multi- 
plicity (exclusive in the 'smeared' sense discussed above, that is, exclusive above scales Qhad = l/^had 
but inclusive for smaller scales and hence still infrared safe). The first underbraced line describes the 
shower-improved contribution to the distribution of O from 'events' which have no perturbatively re- 
solvable emissions at all, the second contributions from events which have exactly one resolved emis- 
sion, and so on. 

This expression can now be expanded further, to any fixed order in the coupling. Each A contains 
one power of Os, and the exponentials inside A must also be expanded. The latter expansion gives rise 
to higher-order coiTcctions which do not increase the parton multiplicity, and thus correspond to the 
'virtual corrections' generated by the shower. The expUcit A factors, in contrast, represent the shower 
approximation of corrections due to real radiation. 



There are several possible definitions of what one might mean by 'matching', reflecting the general 
concept of making two different asymptotic expansions of a given observable agree in an intermediate 
region; our first task is thus to establish a clear nomenclature. In a perturbative calculation, observables 
will have expansions in the strong coupling a^. Each observable will start at some order j in the cou- 
pling, and suffer corrections at subsequent orders. As is usual in fixed-order calculations, we will refer 
to a calculation of the first order for a given observable as being of leading order (LO), a calculation 
accurate to the following order as next-to-leading order (NLO), and so on. Because of the presence of 
infrared singularities, a given order in an observable will receive contributions from perturbative ampli- 
tudes (matrix elements) of different loop order. We will label perturbative amphtudes of 'bare' partons, 
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complete with their infrared singularities, by loop ordei[H That is, we reserve the nomenclature LO, 
NLO, etc. for observables only, and that of tree-level, one-loop, etc. for matrix elements only. 

Event samples as produced by a matched Markov-chain parton shower can be used to measure many 
different observables. We therefore believe it would be misleading to chai^acterize them as being of 
leading or next-to-leading order. More properly, they should be characterized in terms of which matrix 
elements are included in the matching. In this paper, for an arbitrary shower initiator process X (the 
parent process), we intend that "X matched to X + n parions at tree level and X + m parions at loop 
level" (with m < n) should fulfill the following: 

• It should resum the leading soft and coUinear logarithms to all orders, that is, it should be accurate 
up to subleading logarithms. 

• For any j < n, it should reproduce the LO distribution of any observable whose expansion starts 
at order a^, up to connections of order ai~^^ and/or Qhad/Qx' where Qhad is a hadronization scale 
and Qx is a hard scale associated with the X process. 

• For A; < m, it should also reproduce any such distribution calculated at NLO, that is up to correc- 
tions ofO(«^2)and QLdM- 

The first point corresponds to the pure parton shower, the second to CKKW- or tree-level matching, 
and the third to a generalized variant of loop-level matching of the sort implemented by Mc@NLO. 
Our purpose here is to combine all three into a unified approach, in which essentially any tree-level or 
one-loop matrix element could be incorporated with a minimum of effort. 

As an example, consider Z decay. A pure parton shower resums the leading logarithms, but will only 
be "matched" to Z + partons at tree level, as it will generically introduce eiTors of 0{as) in any ob- 
servable. If we match to Z + 1 parton at tree level, we will correctly reproduce any large-logarithm-free 
distribution as predicted using the Z qgq matrix element (the only three-parton leading-order matrix 
element). If we match to Z + 2 partons at tree level, we will correctly reproduce four-jet distributions 
with only ©(af) corrections, con^esponding to the use of tree-level four-parton matrix elements such as 
Z — > qggq. Tree-level matching up to three additional partons combined with one-loop matching up to 
two additional partons, will allow us to reproduce four-jet distributions up to corrections of 0(0^), in 
addition to resumming the leading logarithms, and so forth. 

Let us fix an (ai^bitrary) observable O as representative of the distribution above, and consider a 
computation of the cross section, or partial width as the case may be, differentially in O. We seek a 
prescription that will yield a generated event sample from which distributions can be made that simulta- 
neously fulfill all the three requirements above. 

To specify our matching prescription, we introduce two finite matching terms at each order in Ug, one 
for resolved radiation (R) and one for single unresolved and one-loop corrections (V). Our expression 
for the matched-shower-improved (MS) distribution is then: 



da 



MS 

k=0 



n „ 

x+k)S{{p}x+k,0) , (17) 



where w^x+k tree-level matching coefficient for X + k partons and w^^^j^ is the corresponding 

virtual one. Denoting the couplings present in the parent process X collectively by ax, these matching 



^^We will, however, leave the use of and details of dimensional regularization and infrared cancellations implicit. 
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terms are of order axd^ and axc^s^^, respectively. One or both of them may be zero in the absence of 
matching. The operator 5 is the same as above, embodying an all-orders resummation of both real and 
virtual corrections in the leading-logarithmic approximation. 

The tree- and loop-level matching terms w^-^]^^. and may now be derived by expanding the 

real and virtual terms of 5 separately and, order by order, comparing the contribution from each fixed 
parton multiplicity to the observable O as calculated in a fixed-order expansion: 



da 



ME 



k=0 



X+k 



6(0 - Oi{p}x+k) 



(18) 



where k still represents the number of legs and i represents the number of loops. It goes without saying 
that matching at incomplete orders will involve some arbitrariness, which will be explored further below. 



3.1 "Matching" to X + partons at tree level 

Matching to "X + partons at tree level" just means verifying that the lowest-order expansion of the 
shower is identical to the lowest-order parent matrix element, i.e., that all corrections generated by the 
shower are of higher order. This is trivially true, but let us verify it explicitly as a first exercise. Only the 
first line of eq. ([T6l) is relevant, with the Sudakov A expanded to unity. 



da 
dO 



PS 



/ d«5x wfl, 6 (O - Oap}x+o)) , (19) 



from which we infer that the trivial condition ty^+o ~ l-^x+oP ensures that the two descriptions match 
at lowest order. 



3.2 Matching to X + 1 parton at tree level 

At order as, the parton shower (PS) generates only one term that contributes to X + 1 parton (cf. the 
matching definition above): 

PS : J d^x ImJ'IoP £"f ^ / ^fe^ " - ^(M^+i)) 

(20) 

where the preceding (trivial) matching has been used to replace w^x+o l-^^x+ol^ above, i 

runs over the discrete different branching possibilities, and Ai contains the explicit factor a^. To match 
this to the tree-level {X + l)-parton matrix element, we first divide the complete matrix element phase 
space into a 'resolved' part at early (perturbative) times t < thad corresponding to the region populated 
by the shower, and an unresolved part t > thad which we will treat later. In the region with resolved 
perturbative radiation, we will compute the difference in the distribution of O by subtracting the shower 
term from that of the relevant tree-level matrix element (ME), 

ME : / d$x+i |mJ'|j2^(0-0(Mx+i)) . (21) 
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Subtracting the parton shower term, the matching term (MT) at this partial order, differentially in the 
additional parton's phase space, is simply 

\mP+,? - E 0(*'"(M^+i)-*nM^(-)|Mi°|o(te}x)P, (22) 



ME 



PS 



where i now runs over the number of possible contributing 'parent' configurations for the phase space 
point in question and the G function accounts for the shower not producing any jets harder than its 
starting scale (cf. e.g. the discussion of 'power' vs. 'wimpy' showers in ref. [33]). At such scales, the 
shower effectively has a dead zone, and hence the matching term becomes just the unsubtracted matrix 
element. We show this mostly to illustrate the principle. Because the theta function only affects the 
hard emissions, and because the antenna radiation function captures all the infrared singularities of the 
matrix element, the subtracted matrix element is finite in all single-soft or simple coUinear limits. 

For final-state showers, we can start the shower at a nominally infinite resolution, Qh = l/tn ^ oo, 
i.e. at t = 0, obviating the need for the explicit G function (we defer a discussion of initial-state showers 
to future work). We shall assume this to always be the case and thus define the equation for the matching 
term by: 

MT : w^^l = \mP_^,\' - J: Ai...)\MP^,mx+o)\' . (23) 

As before, the hatted momenta pi appearing in Af j^^Q are a shorthand for the momenta obtained by 
operating on the {X + l)-parton configuration with an inverse map of type i, 

{p^}x = {'^;\{p}x+l)}x■ (24) 

From this discussion, it becomes clear how important it is that the phase space map allows a rel- 
atively clean phase space factorization such that the nested sums and integrals in eq. (l20l ) produce a 
manageable number of subtraction terms with simple borders. Below, we will construct the ViNCiA 
showers explicitly with this goal in mind. 

Before considering the unresolved and virtual corrections, let us remark on a few noteworthy as- 
pects which appear already at this level. We noted above that leading-log resummation only fixes the 
soft/collinear singular terms of A, so that variations in its finite terms are a source of uncertainty for the 
shower, and indeed can be used to estimate these uncertainties. We see here how the matching explicitly 
cancels such variations and hence reduces the uncertainty: if A is made "harder", then the shower gen- 
erates more branchings, but the subtraction term in the matching equation also becomes lai^ger, making 
the matching term smaller and compensating the change. 

An extreme case arises if A is made so large that the matching term becomes negative in some region 
of phase space. This just means that the shower is overpopulating that region relative to the matrix 
element, and hence a negative coiTcction is needed to counter-balance it. The corresponding correction 
events would have negative weights, but there is otherwise nothing abnormal about such a situation. 
Alternatively, one could switch to a shower re-weighting procedure as done in Pythia and thereby 
maintain positive event weights, but in the interest of simplicity we shall not consider reweighting in 
this paper. 

The overall normalization of the parent process under study does change, by an amount given by the 
integral over the matching term. This includes an integral over the arbitrary finite terms in A. However, 
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as we have yet to fix the corresponding virtual term, this just coiTesponds to changes within the range of 
tree-level uncertainty. 

More importantly, this subtraction should in principle be easy to automate. Given any tree-level 
matrix element, which these days can be easily obtained from standard tools like CompHep/CalcHep 
[34], MadGraph [19], and others, the only additional ingredient needed is a subtraction term, whose 
most general form is a sum over lower-point matrix elements multiplied by evolution kernels. As men- 
tioned above, because the leading singulaiities of the resulting subtracted matrix element are absent, it 
should be substantially easier to integrate efficiently over phase space than its unsubtracted counterparts. 

Finally, we note that the matching scheme described above is inherently incremental. With presently 
available methods, a sample of unmatched events cannot easily be modified to produce a matched sample 
(except by doing sophisticated reweightings). Instead, a complete new sample must be generated using 
the matched generator. With our method, a pre-generated set of events need not be re-generated to 
improve the matching; we need only generate an additional set of events corresponding to the matching 
term in eq. (1231 ) and add it to the first, with the relative weight of the two samples fixed by the relative 
integrated cross sections and the number of events in each sample. Of course, this only works to the 
extent that the particular A chosen for the first sample is known at the time the second set is generated. 
The total cross section corresponding to the combined sample would again be different than that of the 
original parton-shower sample, but the difference is of higher order. 

3.3 Matching to X + partons at loop level 

To include the full 0{as) coiTcctions to the initiator process X, and thereby fix the normalization of 
inclusive {X + 0)-parton observables (such as the inclusive cross section or the total width) to NLO, we 
now turn to the relative order as corrections to X with zero additional (perturbatively resolved) partons. 
Again, the parton shower only generates one term, the first term in eq. (fT6l ) with the Sudakov (eqs. ^ 
& (fTTT) ) expanded to order Os : 

H] 

d^x \mP\H{o - o{{p}x) lyx+iYl J - t^H{p}x+i))M...) 

(25) 

which, as a consequence of the unitary construction of the shower, is essentially identical to the real 
radiation term in eq. (l20l) . The only differences are an overall minus sign, and the fact that the observable 
is here evaluated on the parent configuration (X) rather than on one with an additional emission (X + 1). 

We now wish to find the matching term that, together with the one for real radiation above, will give 
the full 0{as) coiTcctions, possibly modulo power coiTcctions in the non-perturbative cutoff thad- To 
accomplish this, we need to include two terms from fixed-order matrix elements, one corresponding to 
genuine one-loop corrections and another corresponding to the real radiation below the hadronization 
cutoff, which was left out above: 

ME : / d^x+i\MP^^\^6{0-0{{p}x+i))+ f d<^x2Re[MPM^^^*]6{0-0{{p}x)) . 

(26) 

Let us re-emphasize that the extra parton in the first term is here unresolved (inclusively summed over) 
and hence the observable cannot really depend on it, up to an overall power correction. Within the 
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required precision, the observable dependence is thus the same for all terms in eqs. (|25] ) & (l26l ). which 
we use to justify lumping them together below. 

The matching term will again be defined by the remainder when subtracting off the parton shower 
contribution from the full matrix element. Differentially in d^x the matching term becomes: 

MT: = 2Re[Mf Mi^^)] + |Mf I^W ^M-) + [ dW-?i, 

^Jallt <i<Px 

where we again used the properties of a clean phase space factorization and extended the definition of the 
subtracted wx+i from eq. (|23] ) into the unresolved region. Because the matched matrix element is free 
of soft or coUinear singularities, the last term is just a power correction, below our required precision. 

Note that a © function restricting the shower term to contribute only after tx has again been avoided 
by letting the shower populate the entire phase space. In one-loop matching with additional partons in 
the final state, a theta function similar to that in eq. (OTI ) will be present. We defer a detailed discussion 
to future work. 

As usual, the first two terms in eq. (|27T ) are separately divergent and a regularization must be intro- 
duced before their (finite) sum can be evaluated. The divergences, which are universal, are usually reg- 
ulated using dimensional regularization and the cancellation can be performed in a process-independent 
way. Only the finite terms must be computed anew for each new process. We thus believe that this part 
could also be automated fairly easily, once the required one-loop matrix elements become available. 

We see here how the NLO normalization of inclusive observables is fixed. In the matching of the real 
radiation term above, the LO normalization changed by the integral of w^^\ a quantity which depends 
explicitly on the finite terms in Ai. The same variation is subtracted in eq. (|25] ). The final normalization 
should accordingly be stable up to higher-order and non-perturbative power corrections. 



3.4 Matching to X + 2 partons (and beyond) at tree level 

So far, we have discussed a matching prescription similar to that of the program Mc@NLO, though we 
have here attempted to develop a formalism more readily applicable to the treatment of non-collinear 
ambiguities and associated uncertainties. The next step in reducing these is to include further informa- 
tion from tree-level coefficients deeper in the perturbative series. 

As mentioned in section |2l we shall now limit our attention to evolution variables which fulfill the 
Mai^kov property in the strictest sense, i.e. which do not have any explicit memory of the event history. 
It then becomes irrelevant whether a particular {X + 1) -parton configuration was obtained by parton 
showering from X + partons or from the tree-level {X + 1) -parton matching term. With a uniquely 
defined "restart scale" tx+i in both cases, the subsequent evolution also becomes the same. 

In fact, the Mai^kov property solves nearly the entire problem for us. We are interested in the relative 
order double real radiation term {X + 2 partons) from a shower which we assume has already been 
matched to X + 1 partons above. By virtue of this prior matching, the total {X + 2)-parton contribution, 
for a history-independent evolution variable, is just given by the parton shower off the tree-level {X + 1)- 
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parton matrix element, here differentially in {X + l)-parton phase space: 

PS : \mP_^,\^ f"%x+2Y.^r^^itx+2-t{{p}x+2))Ai{.^^^^ , (28) 

Jtx+i i a*x+i 

which, apart from the replacement X — > X + 1 and the restriction that t not depend on i (the Markov 
property), is identical to the expression in eq. (l20l) . The second-level matching term would then be. 



l^xW = \mP^2\' - T.M■■■)\MP+l\'Q{t^H{P^}x+2) - tx+i) , (29) 

i 

where the function expresses that the shower evolution is ordered, i.e. that tx+2 must come after 
tx+i- That is, the matching coefficient w is obtained precisely by subtracting the leading singularities 
as expressed by the evolution kernel, along with whatever finite terms we have chosen to include. It 
is essentially the same as the subtraction term (up to the hard-emission modifications due to the theta 
function) that would be used for real emission in a next-to-leading order calculation. 

Note that the correction term for configurations which cannot be obtained from a sequence of two 
ordered branchings is thus unsubtracted. This is most obviously the case for subleading color topologies 
that can arise at higher orders, like Z — > qggq with the gg pair in a color singlet state, but in more 
generality the correction term is simply the full matrix element for any configuration for which the 
parton shower is zero. 

Because of the Markov property, this procedure can be repeated for tree-level matrix elements with 
an arbitrary number of additional emissions. However, while the LL antenna functions still only contain 
the leading singularities, the full higher-order matrix elements will generally contain sub-leading singu- 
larities as well. This leads to problems with unwanted contributions coming from matrix elements with 
"too many" final-state partons. In general, all the following terms may appear (after integrating over 
phase space), 

2n 

< n inmln'"-"^(Q?/Qi), (30) 

m=0 

where Qi 2 are scales in the problem and Lm„ are finite coefficients. For example, at each order n, 
m = is the double logarithmic (eikonal) term and m = 2n is the non-logarithmic (" finite") one. 
We can now be more specific. Since the shower generates only the (leading-color) m = and m = 1 
pieces exactly, the subtraction in eq. (l29l ). beyond n = 1, may leave pieces inside the matching term 
which would be divergent were it not for the hadronization cutoff. If left alone, this would lead to 
distributions of physical quantities with overly large subleading log contributions (divergent in the limit 
the hadronization cutoff is removed), which is obviously not desirable. 

At n = 2, coiTcsponding to {X + 2)-parton matching at tree level, these divergences would be re- 
moved in an NLL shower (where "next-to-leading" here means with respect to the LL shower). Though 
we do make some remarks aimed in this direction at the end of the paper, we note that even if we were 
able to present a complete solution, the same problem would then just appear at NNLL level when at- 
tempting {X + 3)-parton matching, and so on. To do tree-level matching beyond one additional parton, 
clearly, we need a prescription to consistently regulate the subleading logarithms in tree-level {X + n)- 
parton matrix elements. The uncertainty they induce is within the stated accuracy of the calculation, as 
they are higher order in both logarithms and powers of the strong coupling. 
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One possibility is to nominally subtract the subleading logarithms as well in eq. ( |29l ). to the ex- 
tent they are known. Although the LL shower wouldn't regenerate them, this procedure would at least 
cure the problem without affecting the validity of the approach, up to subleading logarithmic correc- 
tions. However, for tree-level matching to X + many paitons, the analytic form of all the corresponding 
jyjmanyj^]^ tcrms would then have to be explicitly subtracted, clearly overkill considering that all we ai^e 
really after is just a regulator. 

A simpler approach is to place explicit restrictions on the 'Wx+2 phase space, cutting out the regions 
where the subleading logarithms become important, for instance by introducing cuts on parton-parton 
invariant masses or transverse momenta. As a rule of thumb one should probably choose the cut to be 
much smaller than the hard scale Qx (so as not to disturb the matching in the hard region) but still 
sufficiently large that ln(Q^/(5^^^) is not much greater than unity. A back-of-the-envelope estimate 
would be that roughly one order of magnitude between the two scales could be a reasonable starting 
point. 

Finally, an alternative approach is to use a Sudakov or Sudakov-like function as a regulating fac- 
tor. This smoothly suppresses unwanted configurations while simultaneously maintaining a fixed-order 
expansion that begins at unity over all of phase space. 

This gives the following general form for leading-order matching with any number of additional 
paitons 

MT : wfl^^ = A({p}x+i+„)(|M(°) 



X+l+n - ^K\l''IX+l+n)y\^y^x+l+n 

-5]A,(...)|Mj°|j2e(tH(te}x+n+l) -tx+n)) ; t < thad , ^^^^ 
i 

where A is either the Sudakov-like function just mentioned or, alternatively, just a G function for the 
cut-off case mentioned above. 

For automated approaches, the 6 function method is probably more appropriate for a stand-alone 
matrix element generator, which would not have the shower Sudakov readily available, whereas a more 
integrated solution could more easily make use of the smoother Sudakov suppression. 

How does this work in practice? To generate a sample of events matched to n additional partons at 
tree level, we should generate events with zero through n partons according to probabilities given by the 
subtracted matrix elements of eq. (|3TI ). and then evolve each configuration using the parton shower. 



4 Non-perturbative Corrections 

The traditional approach in Monte Carlo parton-shower generators is to cut the shower evolution off at a 
low value of the evolution scale, Qhad of C(l GeV). At this point a transition to a different "evolution" 
is made, in the form of QCD-inspired phenomenological hadronization models which explicitly enforce 
confinement and other non-perturbative features. 

From the point of view of a perturbative calculation, this cutoff is simply an arbitrary infrared regu- 
lator, below which partons are not resolved. In the context of the ordered evolution of parton showers, 
however, it represents a scale at which non-perturbative components of the evolution become significant, 
and hence at which point the perturbative evolution kernels used in the parton-shower approximation no 
longer suffice to describe the physics of events; that is, the "evolution" should really contain large cor- 
rections e.g. from pion resonances. 
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In the context of the ViNCiA code, we start by defining the infrared cutoff Qhad in a more uni- 
versal way. Because it simply represents a separation between regions with "large" and "small" non- 
perturbative corrections, respectively, it is not necessary to tie it to the perturbative evolution variable. 
Any infrared-safe phase space contour will do. For instance, one could easily imagine defining a 
hadronization cutoff in terms of dipole-antenna masses applied to a shower which uses transverse mo- 
mentum as its evolution variable, as long as the former regulates all perturbative divergences and the 
latter separates off all regions where hadronization corrections are expected to be large. 

We denote the phase space contour defining the hadronization cutoff for an n-parton configuration 
by 6had(M 

^ /riN ^, /riNN |l in "perturbative" region 

Q...{{p}n) = en,., - t,.,i{pU)) = I Q ^^Ln-perturbat J' region ^'^^ 

where thad = 1/ Qhad is the value of the hadronization cut-off and its functional form (which may be 
different from that of the evolution variable) is given by thadCip})- The Sudakov factor then takes the 
form, 

nf ftcnd r d^>|^L 
exp - / dt„+i / — S{tn+l - t({p}n+l))©had(Mn.+l)^i(-" 

(33) 

for brevity, we have rewritten the sum over i in eq. Q in product form. The perturbative shower 
termination scale tend can now be taken to infinity without any problem, as the divergences are explicitly 
regulated by 0had- The probability that the configuration emits no perturbative (resolved) radiation at 
all is, 

- / dtn+l / -^^^6{tn+l - t{{p}n+l))Qha.d{{p}n+l)Ai{...) 

(34) 

corresponding to Qcnd = l/^cnd 0. This probability is non- vanishing. The matching equations in 
section [3]remain unaltered by the introduction of this hadronization cut-off, except for the replacements 
*end ^ cxD in the integral boundaries along with A — > ©had^- 

The hadronization cut-off has traditionally been imposed in terms of the evolution variable itself, 
since, getting one job done well, it usually gets the other done almost as well. (A few additional cutoffs 
are normally imposed, e.g. to avoid systems with very low invariant masses, but those are minor points). 
This has the disadvantage of making the region defined to be "non-perturbative" different from shower 
model to shower model, and hence a hadronization model fitted with one shower cannot be used as is 
with any other shower 

Decoupling the form of the hadronization cut-off, as proposed here (l34l ). from the shower parameters 
(and in particular the evolution variable), would make the non-perturbative modeling more universally 
applicable. This should be true up to the uncertainty inherent in the perturbative evolution itself. 

This would also be a step towards making it meaningful to compare different parton showers before 
hadronization. This is in stark contrast to the present situation, where different paiton showers are far 
from directly comparable, each having its own cut-off along its own contour Fixed-order parton-level 
calculations could then be replaced by parton showers not including hadronization and matched to fixed 
order matrix elements as the "gold standard" of what is a good perturbative QCD calculation. 
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Figure 1: Left: illustration of the two original dipole antennae in a closed color-singlet gg system in 
the center-of-mass. Right: the system after one branching, showing the branching phase space invari- 
ants Sar and Srb- The (p angle corresponds to rotations around the axis of the original dipole. The ip 
angle corresponds to a rotation of the branching system about an axis perpendicular to the branching 
plane. The "forbidden angle" 9, always set to zero in our maps below, would correspond to rotating the 
branching plane off axis with respect to the original dipole. 



5 The ViNCiA Code 



We now turn to a proof-of-concept implementation of the ideas contained in previous sections, in the 
form of the ViNCiA code (Virtual Numerical Collider with Interleaved Antennae), implemented both as 
a stand-alone program and as a final-state shower plug-in for PythiaS. H — > gg, matched to i7 — > ggg 
at tree level and H ^ gg a.t loop level, has also been implemented in both versions, according to the 
matching terms defined in section |3] For the plug-in, this includes the possibility of generating negative- 
weight connection events when the shower is overpopulating phase space. 

The numerically implemented shower is based on an interleaved evolution (see e.g. ref. [3]) of 
systems of color-ordered QCD antennae. The implementation discussed here is limited to gluons, and 
uses a strict dipole-antenna factorization [28]. (The name 'dipole factorization' is associated with a 
related NLO formalism due to Catani and Seymour [35].) Inserting the massless dipole phase space, 
eq. (fTOl ). into the event Sudakov eq. (l34l) yields a product of individual color-ordered dipole Sudakov 
factors 



A(tn, tend', {p}n,thad) = 

[[exp i- J^dSar ^^rb ^(5(^1 " *(Mn+l ))6had (Mn+1 ) J^^j^^ j , 

(35) 

where the branching invariants Sar = {Pa + Pr)^ and Srb = {Pr + Pb)'^ are illustrated in Fig.H] We now 
proceed to give explicit forms for each of the objects required by section |2]for the definition of a shower. 



5.1 Evolution Variable 

We shall here consider only Lorentz-invariant evolution variables, Qe = We have implemented 
two different choices, corresponding to ordering in transverse momentum and in dipole mass (~ parton 
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ya,=Sa/S yar=Sa/S ya,=Sa,/S 

Figure 2: Illustration of the progression of the evolution variables over phase space for massless partons: 
contours corresponding to yi n =0.5(left), O.l(middle), and O.Ol(right). On the x and y axes respectively 
are the rescaled invariants y^j. = Sar 

/ s and y^b = Srb/s after the branching. The soft singularity is at 
the lower left corner, and the coUinear singularities lie along the axes. 



virtuality) respectively: 

type, . 0, . 4— =4p,^,„„,, ^^^^ 

type II : Qjj = 2mm{sar, Srb) 

where the normalizations have been chosen so that the maximum value of the evolution variable is 
always the dipole-antenna invariant mass s (to avoid cluttering the notation, we now let the superscript 
[i] be implicit). We will usually work with dimensionless versions of these invariants. 




2 Ql A '^ar^rb , 
yf = = 4 2— = ^VarVrb 



yE = { J, , (37) 

: Vn = — = 2min(yar,yrb) 



s 



where y-i = Si/s, so that the maximal value of yi u inside the physical phase space is unity. A com- 
parison of iso-y contours for these two variables in the branching phase space is shown in Fig.|2l Their 
complementary nature is now readily apparent. The transverse-momentum or Ariadne variable (type 
I) will categorize a hard but collinear branching (close to one of the axes) as harder than a wide-angle but 
soft one (close to the origin), whereas the the virtuality ordering (type II) will tend to do the opposite. 
This affects which regions act to Sudakov suppress which branchings during the evolution. 

Note, however, that the definitions in eq. (|37] ) do not yet completely obey the Markov condition. Be- 
cause gluons are indistinguishable, it is not possible to single out the radiated parton r without knowing 
the branching history of the configuration. In other words, when showering off a three-gluon configu- 
ration with an unspecified history (e.g. from the three-gluon matching term), we have several possible 
choices of what "restart" scale to choose, depending on which of the partons we decide to call r. We em- 
phasize that this is not a problem for matching to first order (X + 1 at tree level and X + Oa.t loop level), 
since the history-dependence only has to do with what restart scale to choose and hence, at the earliest, 
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affects the second emission. There is no fundamental difficulty in defining variables which strictly obey 
the Markov condition, but as already discussed in section [3l we postpone a detailed discussion of this 
aspect to future work. 



5.2 Phase Space Map 

We must next choose a phase space map. The restriction to dipole-antenna phase space factorization, 
eqs. & ([Tol l, already fixes most of the 5 functions: all partons except the two involved in the splitting 
are just "copied" to the (n + 1) configuration. Denoting the branching antenna pair by [i], eq. Q implies 



Pj = Pj V j i [i] . (38) 

The branching antenna pair, denoted o and 6, are replaced as shown in fig. [B by a trio of partons, denoted 
a, 6, and r. This replacement conserves energy and momentum, and keeps all partons at their physical 
masses. In the center of mass frame of the parent dipole, the energies are related to the branching 
invariants as follows. 



2^ 



Eb = ' , (39) 



2v^ 



Our discussion here will focus on massless partons, nia = mi, = nij. = 0. Since the phase space 
consti'uction implicitly uses 5 functions requiring the partons to be on shell, the absolute values of the 
momenta are equal to the energies for massless particles. 

Staying in the dipole's center-of-mass frame (DCM), there is no freedom left to choose the relative 
angles between the three daughter partons in the branching plane: 

2EaEr + ml + ml- Sar 

COS Oar = ^1 n \ ' 

2|Pa||Pr| 

2EaEb + nil + ml- Sgb 

cos Oab = 77. n i • (40) 

2|Pa||Pfe| 



We still need to fix the orientation of the three daughter partons with respect to the parent dipole. 
This involves three Euler angles. As noted in section |2j one of these is fixed by requiring that the 
branching plane contains the dipole axis, or equivalently that the normal to the plane be orthogonal 
to the dipole axis. (This imposes parity conservation on the 2^3 transition.) Another angle is just 
the integration variable representing rotations around the dipole axis. While the latter is here chosen 
isotropically we note that the matching terms will still (re-)introduce anisotropics up to the order of the 
matching. 

The remaining ambiguity in the phase space map thus rests entirely with the last Euler angle, the 
one corresponding to rotations ai^ound an axis perpendicular to the branching plane. In the context of 
Ariadne [17], a choice was made which "least disturbed" neighbouring dipole^. In order to explore 



'*For quark antennae, Kleiss has shown that an optimal choice exists [32], but for gluon antennae the situation is less clear. 
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the consequences of this ambiguity, we have so far implemented three discrete possibilities for this 
angle, defined by eq. ^ to be the angle between parton a and the original parton a in the DCM frame, 

V'Ariadne = El + El ^"" ~ ^"''^ ^"^^^ 
^ps = 1 ° ' ^ ] (42) 

[ - t)ab ; Sar < Srb J 
V'Ant = 1 + -^ , (43) 



1 - Vrb 



where yaa in the last Une is defined by: 



/ = , (44) 

Var I Vrb 



p = \ll + ^f{l- f)yaryrb/yab , (45) 

(1 - P)yab + 'ifVarVrb 

Vaa = TTT, ^ • (46) 

2(1 - Var) 

It is important on physical grounds that — > smoothly as parton r becomes colUnear with parton 
b (that is, as Srb — > 0), so that parton a becomes aligned with parent parton a, and likewise when the 
roles of a and b and a and b are simultaneously interchanged. This ensures that the daughter system 
approaches the parent one in this limit. Otherwise, however, there are no constraints on ip. All three 
alternatives satisfy this constraint. 

The first choice corresponds to the ARIADNE map, just as our first evolution variable corresponds 
to the Ariadne one [17]. We hope this helps make comparisons between the two approaches simpler. 
The second corresponds roughly to conventional parton showers, in which the non-radiating parton only 
recoils longitudinally. (Since our antenna shower does not maintain a clear distinction between which 
parton radiates, the one with the largest invariant mass with respect to r is chosen to play the part of 
recoiler.) The last choice is an example of a more general form; different choices of / could be used to 
explore it more fully. 

5.3 Radiation Function 

We have thus far implemented only the gg — > ggg radiation function, for which we have used the 
Gehrmann-Gehrmann-De Ridder-Glover "global" antenna function/3 [36iE 



f3iPa,Pr,Pb) = ^ 



(1 - Var - Vrb) { \ \ I + ^ 

VarUrb Vrb Uar J "J 

"soft" "colhnear" 



(47) 



In this formula, s'*] is the mass squared of the dipole-antenna. "Global" means that the phase space of 
each antenna is unrestricted by overlap with other antennae; the normalization and singularities are such 
that the sum of contributions has the desired structure. The leading (double logarithmic) singularities 



^ Note that we have changed the non-singular term from 2/3 to 8/3, relative to the original paper. 
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correspond to two invariants vanishing (soft radiation), and arise only from one antenna. The single- 
logarithmic (coUinear) singularities receive contributions from two neighboring antenna. 

We choose a second-order polynomial in the invariants for the form of the arbitrary finite terms, 
imposing only the restriction that the antenna function be positive definite. Combining /g above with 



the normalization implied by eq. (1351 ). the radiation function for the ViNCiA gluon shower becomes: 



MPa,Pr,Pb) 



Vrb) 



UarVrb Vrb Var 



a,l3>0 



(48) 



where finite terms are parametrized by the constants Ca/s- We can explore systematically the conse- 
quences of making the radiation function harder or softer by varying Cq,^; e.g., the special case coitc- 
sponding to the antenna function can be obtained by choosing Cqo = 8/3. As discussed in section 
121 matching absorbs these variations in the matching terms, leaving only the uncertainty due to genuine 
higher-order terms in the shower We can thereby quantify the reduction in the associated uncertainty. 

5.4 Renormalization Scale 

We let the renormalization scale for at each branching be given either by the evolution scale at the 
branching, Q^, or by the invariant mass of the dipole being evolved, s^. 




An 



KrQe 



(49) 



where Qe the evolution variable and we allow for an arbitrary prefactor Kr to be applied. For 
example a factor Kr = 1/2 applied to the type 1 evolution variable would yield a renormalization scale 
equal to the Ariadne definition of transverse momentum. By default we will use a one-loop running 
as, but we leave open the option of studying fixed coupling or two-loop running as welj^ 



5.5 Starting and Ending Scales 

For a parent process producing two partons in a decay, such as H gg which we shall consider below, 
we choose the initial starting scale to be the full phase space, s, so that the shower does not have a 
dead zone. After branching at scale Qe, the shower evolution continues from that scale. As already 
discussed, this does imply a slight dependence on the shower history, as the same configuration can in 
principle be obtained by different branchings coiTcsponding to different values of Qe- For showers off 
the thiee-gluon matching term, which has no history to provide a unique scale, we compute the scale 
corresponding to each possible ordering and select the smallest of these, as the matching is intended to 
describe the hardest emission. The shower is cut off in the infrared by an evolution-independent contour, 
as described in section IH The choices possible for the functional form of this contour are the same as 
for the evolution variable, eq. (l36l ). The history dependence could be eliminated by using an antenna 
function restricted to a 'wedge' or 'sector' of phase space; we leave further discussion of this to future 
work. 

*In the PythiaS plug-in, we rely on the as implementation in Pythia, which likewise provides these choices. 
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5.6 Shower Implementation 

Shower generation proceeds as follows. Given a starting scale Q„, a trial branching for each antenna 
dipole is found by generating a random number R G [0, 1] and solving for Qn+i in the following "trial 
equation": 



R 



^{Qn, Qn+l) 

rQn 



exp 



dQn+l / dSa-, 
'n+l JO Jo 



dSrb S{Qn+l - Qi{p}n+l)) 



A{.. 



2Ji] 



(50) 



where we use A to signify that a nominally larger branching probability A> A may be used to generate 
these trials (for instance using an over-estimate of a™^^ > as{fJ-R), no hadronization cutoff, etc); the 
resulting distribution will then be corrected by subsequent vetos. 

In traditional approaches, an equation for Qn+i{Qn, R) is obtained by analytically inverting eq. (l50l ). 
Since we wish to be able to choose arbitrary evolution variables and radiation functions, however, we 
have instead used a more numerical approach. 

For fixed coupling, the Sudakov factor only depends on one quantity, the ratio of the evolution scale 
Qn+i to the starting scale Qn- Re-expressing the Sudakov factor in terms of dimensionless ratios of 
invariants. 



A(ytrial) = A(l,ytrial 

= exp 
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ytrial "'0 JO 



167r2 



(51) 



where Ve = Qe i& the dimensionless evolution scale, as defined by eq. (|37] ). Because the combination 
sWA(...) is independent of sW, cf. eq. (1481 ). this quantity depends only on a single variable, ytriai- 
Accordingly, it is simple to tabulate it during initialization; we do so using a cubic spline, performing 
the integrals inside the exponent either numerically (via two-dimensional adaptive gaussian quadrature) 
or analytically (as a counter-check). 

We may then solve the equation R = A (ytriai) numerically for ytriai using the splined version of the 
Sudakov and standard root finding techniques. These are computationally quite efficient. 

The antenna with the largest trial scale is then selected for further inspection. A <p angle is chosen 
uniformly, and the remaining degeneracy along the iso-y contour (as shown for example in fig. ^ is 
lifted by choosing a complementary invariant, which we call z, according to the probability distribution. 
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(52) 



where \ J{yE,z)\ is the Jacobian arising from translating {yar,yrb} to {yE,z} and z^amiyE) is the 
smallest value z attains inside the physical phase space for a given yE- Since z merely serves as a 
parametrization of phase space along an iso-y^; contour, its definition is arbitrary, so long as it is linearly 
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independent of ije- Depending on the type of evolution variable, we choose z as 

type I : . = , ^^^^ 

type II : z = max{yar,yrb) , 

leading to the Jacobian factors | Ji| = l/(4z) and | Jn| = 1/2, respectively, and the phase space bound- 
aries 

type I : ZmmiVE) = 5(1 - ^1 - Ve) , Z^a.yi{yE) = 1(1 + " Ve) , 



type II : z„,m{yE) = hvE , z^nxiVE) = l-\yE ^^^^ 



where the type II case should be divided into two branches, one with yar > yrb and one with y^t > yar, 
each having the phase space limits given here. 

Because the Iz functions depend on two independent variables, yE and z, we have not implemented 
a splined approach for this task. Instead, we use analytical integrals over the two kinds of phase space 
regions we are interested in. In generic form, these are 

max AT 

hiVE, z) = -^—^ [S{z) - S{zr^^) + K{z) - K{zr^^)] , (55) 
47r 

with S{z) coming from the soft and collineai^ singular terms, 

yl yE ^ 8 + yE z ^ -12z + 3z^ - 2z^ + I2ln z 

''''' ■ ^'^'^ = 1^-3^ + ^-4 + 6^E ' 

type II : Su{z) = ^ (-6(8 + y| + 2z - - Sz^ + 3(2 - yi=;)(8 + y|) Inz 

and K{z) coming from the finite polynomial. 
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type I : K,{z) = W^^"" 
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(57) 



Vrb) 



4^V4 
1^„ z^+' 

a,P 

with a,(3 >0. 

With the trial resolution scale yE and the energy-sharing fraction z now in hand, we can compute 
yar and y^b- Together with (chosen above), this gives the complete branching kinematics. We now 
apply a veto, accepting the trial branching with probabiUty, 

D n> /! \ \ n ^ Q^(^fi(M"+i)) .-os 

-Taccept = WhadiVhadUP/n+lj ~ Vhadj • (Jo) 

That is, the branching is only accepted if it is inside the perturbative region and then only with probability 
Us/a^^^, which reduces the effective coupling to the correct value by virtue of the veto algorithm. (Note 
that the event is not thrown away, it is merely the branching which is vetoed.) In order to evolve the 
system further, we repeat the steps above. The trial branching scale becomes the new starting scale, 
whether the above branching was vetoed or not. 

The evolution continues until there is no perturbative evolution space left (the equivalent of reaching 
the hadronization cutoff in our terminology). In the cunent implementation, we consider this condition 
satisfied for a given antenna if ten consecutive trials are rejected due to the 6had condition. 
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5.7 Matching Implementation 



In this paper, we restrict ourselves to matching at first order (at tree and loop level) for a scalar decaying 
into two gluons via an effective point coupling. By first order matching we mean that, in addition to 



W2 



\M. 



(0)|2 



2 I , we include the matching coefficients w\ and W2 ■ For the decay process H ^ gg the 
subtracted matrix elements are relatively easy to obtain. Given the Born squared matrix element [Mg*^^ p 
for H ^ gg we find for the 2-gluon matching term at one loop 
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and for the 3-gluon matching term 
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(60) 



where 



Farb — X! ^a/SUary. 



rb 



o,/3>0 



(61) 



We note that by taking Cqq = 8/3 and all other coefficients equal to zero, the three-gluon matching term 
is zero. This means that the complete H — > ggg amplitude has been absorbed into the Sudakov factor. 



5.8 Preliminary Results 

We now turn to a comparison between results obtained using a few different parameter and variable 
choices with the ViNCiA code (using the ViNCiA plug-in module with PYTHIA8.086). Recall that we 
are here studying pure gluon evolution in the fictitious decay of a scalar to two gluons. We thus intend 
these results mostly for illustration of the method. We use the type I evolution variable with a one-loop 
running Og- The hadronization scale is chosen to be Qi = 1 GeV. 

The plot in fig. [3] illustrates the distribution of the (symmetrized) type-I resolution scale for three- 
parton configurations, obtained with unmatched ViNCiA for "soft" (all Ca/s = 0), "standard" (Cqo = 
2.66 ~ 8/3) and "hard" (Cqo = 10) variants. For all curves, jj,R = Qi/2 ~ p±. The point Q = 
80 GeV corresponds to the "Mercedes" configuration. While the variations greatly affect the shape of 
the distribution, the peak position remains fairly stable, here at around a tenth of the original mass. 

To investigate how matching reduces this uncertainty, fig.|4]shows the two-dimensional phase space 
population for three-parton configurations corresponding to five different settings of the ViNCiA plug- 
in, from top left to bottom right: soft (unmatched), soft (matched), standard, hard (unmatched), hard 
(matched), where matched here refers to matching to H ^ ggg at tree level. As in fig. |2l the dipole 
phase space is represented as a triangle in the two phase space invariants, y-ij = s-ij/s and yjk = sjk/s, 
here symmetrized over i, j, and k (because gluons are indistinguishable). The dark color in the center of 
the plots indicates low probability, with warmer colors (lighter shades) denoting increasing probabiUties 
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Q [GeV] 

Figure 3: Type I evolution scale for the first shower branching, symmetrized over a, r, and b: Q = 
min((5i(a, r, b), Qi{r, b, a), Qi{b, a, r)) ~ 2p±, for H ^ gg with ruH = 120 GeV, showered with three 
different choices of the finite terms in the ViNCiA shower, coiTcsponding to "standard" (Cqo = 2.66), 
"hard" (Cqo = 10), and "soft" (Cqo = 0) variations. 



towards the corners and sides. In order to focus on the hard central region, the color scale has been forced 
to white for l/NdN/dyijdyjk > 2, thus "whiting out" the strong peaking towards the corners and sides 
of the triangle which would otherwise dominate. In the top row, con^esponding to the soft shower, the 
matching fills in missing configurations in the near- "Mercedes" region. In the middle row, no matched 
plot is shown, because setting Cqo = 8/3 corresponds to exponentiating the matrix element itself and 
hence the shower already produces the "correct" result in this case. For the lower two plots, the shower 
is significantly harder than the matrix element. In this case, the code responds by generating a matching 
term with negative weight, effectively reducing the population of the hard region when "added" to the 
unmatched events. The reduction in uncertainty for this observable is evident by comparing the variation 
from top to bottom on the left, with the variation on the right. (The somewhat odd-looking contours are 
merely an artifact of Root's contour algorithm operating on binned histograms.) 

In fig.|5]we show the number of two-jet configurations as a function of the type-I resolution scale, 
roughly equivalent to the Sudakov factor expressed in this variable, before (left) and after (right) match- 
ing to the tree-level H ggg matrix element. As could be expected, the uncertainty on this observable 
is greatly reduced by matching to this level. 

The thrust is a more complicated observable, whose distribution is shown in fig.lH again with (right) 
and without (left) matching. We see that, in the region accessible to 3-parton configurations, 1 — T < 
1/3, the variation is indeed canceled, while in the region of large 1 — T, accessible only to 4-parton 
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Figure 4: Phase space population l/NdN/dyij/di/jk, for three-gluon configurations, symmetrized over 
all combinations of i, j, and A;. Top left: the "soft" ViNCiA shower off two-parton configurations, with 
all Cij = 0. Top right: as in top left, but including the tree-level H ggg matching term. Note how 
the radiation "hole" in the center is filled in slightly. Middle: the "standai^d" ViNCiA shower, which 
absorbs the tree-level H —>■ ggg matrix element correction. Lower left: same as upper left, but with 
Coo = 10, corresponding to a "hard" ViNCiA shower. Lower right: as in the lower left, but including 
tree-level matching. In this case, the H — > ggg correction is negative, reducing the population of the 
central region. 
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k_L[GeV] k_L[GeV] 



Figure 5: The number of exclusive two-jet configurations as a function of Durham k±. Labels are 
identical to those in fig. [3l except that we have here used a different renormalization scale choice, 
I^R = Qii/2, translating to ij-r = m///2 for the first branching over all of phase space. 



and higher configurations, some uncertainty remains (though with 2.5M events generated per curve, 
small statistical fluctuations also become noticeable in the shaiply-dropping tail). We also note that the 
behaviour at small (1 — T) is very sensitive to the choice of fiR, a point to which we plan to return in 
the future. We have used /x/j = Qu/2 here. 

6 Still Deeper? 

In sect. 121 we discussed how to match a leading-logarithmic parton shower to tree-level calculations 
with an arbitrary preselected number of resolved partons. We also outlined how to perform a similar 
matching at one loop. In the previous section, we presented a first implementation of these ideas. While 
there is still a great deal of work to be done in fleshing out and implementing our approach, it is also 
interesting to peer ahead, and ask: how would one go further in perturbation theory? How could one 
further improve the accuracy of parton-shower predictions? 

One can presumably proceed to higher fixed orders, matching to NNLO calculations by deriving 
generalizations of the equations presented in sect. [3l We shall not examine such matching in any greater 
detail. Instead, let us explore the possibility of resumming subleading logarithms, that is including not 
only terms of 0{a^ In^"'^""^ y) where y is a lai-ge ratio of scales, but down to 0(0"+^ jj^2n,2n-i yy 

For this purpose, it is crucial to have a formalism that treats all the leading-logarithmic singulari- 
ties exactly point-by-point in phase space. This is true of the antenna-based formalism described here 
(and would also be true of showering based on the Catani-Seymour dipole formalism). To set up a 
subleading-logarithmic shower, we must consider conections to the showering kernel itself. There ai^e 
two kinds of corrections: virtual corrections, and real-emission ones. The former still correspond to a 
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2 — > 3 branching process, but with the branching probability computed to one order beyond leading 
in as- The real-emission corrections correspond to a new branching process, 2 — > 4 partons. Such a 
branching can occur with the basic Sudakov, of course, but only in two branching steps. Here, it will 
sometimes happen in one step. Indeed, the kernel that will ultimately enter a modified Sudakov factor 
will not be simply the 2 — > 4 branching or antenna function, but rather that function, with the iterated 
2 — > 3 contribution subtracted out. This excess represents the genuine correlated 2^4 branching 
probability. The required ingredients in a dipole-antenna approach — one-loop connections to the basic 
antenna, and the 2^4 tree-level antenna function — are known from the development of an NNLO 
fixed-order formalism [37]. 

We would further need a definition of the evolution variable that can be evaluated on n — > n + 2 
branchings, and that regulates all infrared divergences in them. We would also need an appropriate 
phase-space mapping for the following factorization, 

d$„ d«5„_2 d$2 ' 

which is now six- rather than three-dimensional. With these definitions and mappings, the NLL Sudakov 
would presumably take the form, 

ANLL(in,iond) = ^LL(*n, ^cnd) X 



exp 



/■iond /■ d$'"'! o r 1 r 1 

dtn+2 E y ^|^5(W2-t[^kMn+2))4l4(Wn-Mn+2) 



je{n-»n+2}' 



(63) 

where A'll includes the one-loop connections to the LL kernel. We will face the issue of fixing the finite 
terms in the LL Sudakov factor upon matching to the NLL one, and a related issue of maintaining the 
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positivity of the resulting NLL corrections A^^^. The requirement of maintained positivity is crucial to 
a probabilistic interpretation, such as the Markov chain. 

Assuming these issues can be resolved in a satisfactory manner, the matching prescriptions appear 
to generalize in a straightforward way. At order ag, nothing much changes, except for possible mod- 
ifications to the LL kernels. At relative order a^, the shower now simply produces one extra term, 
corresponding to a direct 2 — > 4 branching. It contains the proper subleading logarithms by construc- 
tion, and accordingly the tree-level matching equation will have the form (differentially in d$x+2)^ 

Wx+2 = |Mx+2|'-E4l3(-)|Mx+lPe(t[^l(Mx+2)-tx+l) 

* r64") 

j 

where ^2^3 includes one-loop corrections and ^2-^4 is the direct 2 — > 4 kernel, as discussed above. 

At least from this perspective, the inclusion of explicit n — > n + 2 branchings poses no fundamental 
problem. In addition, inclusion of matching to one-loop and two-loop matrix elements, which will both 
be modified by subtraction terms, would open the way to event generation at NNLO. 



7 Conclusion and Outlook 



We have presented a new general formalism for parton-shower resummations. The formalism allows 
us to explore both the uncertainties inherent in the parton-shower predictions, and the reductions in 
them possible by matching to fixed-order matrix elements. We keep track of the ambiguities of the 
shower approach away from the soft and coUinear regions, allowing us gain a systematic estimate of 
the associated uncertainties. The quantification of these uncertainties, as well of their reduction by 
matching, is novel. 

We have outlined a general approach for matching to fixed-order matrix elements, based on a sub- 
traction approach which generalizes that of Frixione and Webber. We also presented a specific algo- 
rithm based on antenna factorization and dipole-antenna showers, generalizing that of Gustafson and 
Lonnblad. The formalism is simple and intuitive, but is powerful enough to match fixed-order matrix 
elements at higher multiplicity both at tree- and one-loop level. In this respect, it provides a generaliza- 
tion of both the CKKW and Mc@NLO approaches. The (arbitrary) choice of non-singular terms in the 
shower kernel is explicitly canceled by the matching terms, which allows us to quantify the degree to 
which matching to a given order reduces the uncertainty inherent in parton-shower predictions. 

We presented a generalization of the definition of the hadronization cut-off that would make possi- 
ble a more universal modeling of non-perturbative physics, allowing more meaningful comparisons of 
different parton-shower approaches, as well as the improvement of fixed-order parton-level calculations 
without reference to a specific hadronization model. 

We have developed a proof-of-concept level implementation for matching of gluon showers in the 
decay process H ^ gg including both real and virtual corrections, in the form of the ViNCiA code, and 
have presented illustrative comparisons with and without matching for a few benchmark distributions. 

The next step will be to include quarks and perform a more comprehensive study of both H — > 
gg and Z ^ qq fragmentation, exploring the properties of the ViNCiA algorithm and its relation to 
existing approaches in greater detail. We plan to go into greater detail on various theoretical aspects 
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in a future paper [29]. We believe that it should be straightforward to automate matching for general 
lepton collider and decay processes matrix elements, once the evolution variables are generalized to be 
history-independent. The inclusion of initial-state radiation and matching will be necessary to extend 
the approach to hadron collisions. The formalism outlined here should be sufficiently general to make 
this feasible. Indeed, as we have discussed briefly, we believe it will be sufficiently general to open a 
path to matching and showering deeper into the perturbative regime, both in powers of and orders of 
subleading logarithms. 
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